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ABSTRACT 

We present a new approach to calculate the particle distribution function about rela- 

tivistic shocks including synchrotron losses using the method of lines with an explicit 
finite difference scheme. A steady, continuous, one dimensional plasma flow is consid- 
ered to model thick (modified) shocks, leading to a calculation in throe dimensions 
plus time, the former three being momentum, pitch angle and position. The method 
accurately reproduces the expected power law behaviour in momentum at the shock 
for upstream flow speeds ranging from O.fc to 0.995c (F G (f , fO]). It also reproduces 
approximate analytical results for the synchrotron cutoff shape for a non-rclativistic 
shock, demonstrating that the loss process is accurately represented. The algorithm 
has been implemented as a hybrid OpenMP-MPI parallel algorithm to make efficient 
use of SMP cluster architectures and scales well up to many hundreds of CPUs. 

Key words: acceleration of particles, relativistic processes, shock waves, methods: 
numerical 



1 INTRODUCTION 

Particle acceleration at astropfiysical shocks has been the 
subject of a great deal of research for over 30 years. For 
non-relativistic shocks it has been shown that the diffusion 
approximation can can be used to find accurate solutions 
to the governing partial differential equation(s) in a wide 
range of problems (Drury f983). The diffusion approxima- 
tion arises from the assumption of near-isotropy of the distri- 
bution function in the region around the shock front. How- 
ever, in the case of relativistic shock fronts this assumption 
of ncar-isotropy is no longer vahd and it becomes necessary 
to consider the full angular dependence of the distribution 
function (Kirk & Duffy 1999). This effectively increases the 
dimensionality of the governing equation. By assuming a 
gyrotropic distribution this increase is limited to a single di- 
mension in the form of the particle "direction angle", i.e. the 
angle between the particle's mean (gyro-averaged) velocity 
and the shock normal. For plane parallel shocks the pitch 
angle is the same as the particle's pitch angle. 

Recently, a semi-implicit finite difference, method of 
lines approach was used to examine time dependent par- 
ticle acceleration at a thin (discontinuous), non-relativistic 
shock in the presence of synchrotron and inverse compton 
losses (Vannoni et al. 2009). This work has shown that finite 
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difference methods and the method of lines can be used to 
great effect to solve the equations of diffusive shock acceler- 
ation (DSA) in the diffusion approximation. We would like 
to dispense with the diffusion approximation and the test 
particle approximation but, as we will see, this dramatically 
increases the complexity of the problem and necessitates the 
use of high performance computing (HFC) techniques. 

Perhaps the most successful method that has been used 
in the analysis of DSA at relativistic shock fronts is the 
eigcnfunction method. First published in the late 1980s 
(Kirk & Schneider 1987), the method was initially used to 
investigate the spectral index and angular dependence of 
the distribution function for shocks with upstream Lorentz 
factors of r < 5. Only the steady state solution was con- 
sidered with isotropic, momentum independent diffusion in 
particle direction angle. The effects of shock structure were 
also omitted with a discontinuous shock model being em- 
ployed. By ignoring radiative losses or injection effects and 
assuming the solution to be a power-law in momentum, the 
spatial and angular dependence of the solution were then 
calculated. Soon afterwards Heavens & Drury (1988) inves- 
tigated momentum dependent, anisotropic diffusion by us- 
ing the more general Kolmogorov power spectrum for the 
scattering waves. Following this work, monotonic shock ve- 
locity profiles were considered and the assumption of power 
law solutions was relaxed to allow the investigation of injec- 
tion effects at low energies (Schneider & Kirk 1989; Kirk 
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& Schneider 1989). The assumption of momentum inde- 
pendent, isotropic diffusion was reintroduced for this work. 
About 10 years later an additional paper (Kirk, Guthmann, 
Gallant & Achterberg 2000) extended this method to include 
the full range of shock velocities from non-rclativistic to 
ultra-relativistic. Again, shock structure was neglected and a 
power law distribution in momentum assumed. Anisotropic 
diffusion coefficients, independent of momentum were con- 
sidered, allowing the correlation length of the tangled field 
parallel to the shock to vary with respect to that perpen- 
dicular to the shock. Later high-energy cutoffs due to syn- 
chrotron losses were also incorporated by Dempsey & Duffy 
(2007). 

Given the recent success in the application of the 
method of lines to nori-relativistic shock acceleration, it is 
extended in this work to the relativistic and ultra-relativistic 
cases. Such a method could might possibly be coupled with 
a relativistic (magneto)hydrodynamics code to model mod- 
ified relativistic shocks. For the present work we partially 
retain the test-particle approximation, meaning that we will 
not include the dynamic, time-dependent back-reaction of 
the accelerated charged particle population on the back- 
ground plasma. Wo employ a continuous, steady velocity 
profile describing the shock which can be chosen as an arbi- 
trary monotonic function of position before the calculation. 
Therefore a more "realistic" modified shock structure could 
be imposed but we will not propose such a structure here. 
Moreover, having successfully applied the present methods 
to a steady shock profile, the introduction of a time depen- 
dent shock structure should be relatively straight forward. 
The additional terms required in the transport equation are, 
in fact, remarkably similar those containing the spatial de- 
pendence of the profile and therefore should not produce 
significant numerical challenges. Subsequent coupling with 
a relativistic (magneto)hydrodynamics code could perhaps 
bo considered to study modified shocks. Such a coupling 
would be far from trivial, involving the inclusion of some 
complex physical effects and the associated numerical issues. 
For example, one would have to consider the currents aris- 
ing from the anisotropic distribution of high energy charged 
particles and their effect on the plasma flow. The pressure 
tensor associated with the anisotropic momentum fiux could 
also introduce difficulties. For now we concern ourselves with 
the feasibility of this numerical approach to the acceleration 
problem and we leave these extensions for future work. 

A number of physical challenges have boon presented 
that call into question the validity of Fermi acceleration at 
relativistic shocks. In the ultra-relativistic limit we expect 
shocks to become superluminal with strong cross-field diffu- 
sion becoming a necessary ingredient to support the theory 
(see for example Kirk & Duffy (1999)). The detailed struc- 
ture of the magnetic field and its interaction with the ultra- 
relativistic plasma flow and high-energy particles about such 
shocks is far from clear. Promising results have recently been 
emerging from largo particle in cell simulations which can be 
used to investigate various (micro-)instabilities which can 
scatter particles and significantly amplify or alter the char- 
acter of the magnetic field close to the shock (Spitkovsky 
2008; Dieckmann et al. 2008; Kcshet ct al. 2009). It is clear 
that the magnetic field about relativistic and ultrarelativis- 
tic shocks is quite complex and heavily coupled with the 
behaviour of the high energy particles produced in such re- 



gions. It has been shown that a sufficiently regular magnetic 
field upstream of an oblique, relativistic shock will cause 
particles to leave the loss cone due to the Lorentz force and 
recross the shock having completed < 1 gyration about the 
upstream field lines (Achterberg et al. 2001). In this case, 
pitch angle scattering does not describe the particle mo- 
tion well in the upstream region. Of course the assumption 
of a sufficiently regular magnetic field to support this pro- 
cess may not be valid due to the various instabilities and 
turbulence that may exist in the region. In order to pro- 
ceed with the present method, we will assume that particle 
transport in the upstream region is well described by the 
same equation used to model pitch angle diffusion at par- 
allel shocks. Wo adopt a similar argument to Kirk et al. 
(2000), noting that this equation may well describe parti- 
cle transport in other more complex cases where particles 
interact with a more complicated field. In other words, we 
have assumed that the particles exhibit diffusive behaviour 
in "direction angle" (that between their velocity and the 
shock normal) rather than pitch angle. Therefore, in this 
paper, we consider particle acceleration at plane, parallel, 
relativistic shocks and note that the results may be extensi- 
ble to other more complex and realistic cases. 

We will combine a number of the sub-problems tackled 
separately with the eigenfunction method and seek a numer- 
ical solution. Hence we would like to consider the time de- 
pendent acceleration of particles at shocks with a wide range 
of Lorentz factors (1 < F ^ 10). In addition, we consider 
the effects of shock structure to permit the investigation of 
modified shocks, potentially including time-dependent hy- 
drodynamics. We would also like to include anisotropic, mo- 
mentum dependent diffusion in particle direction angle to 
permit the introduction of a range of physically motivated 
diffusion models. Finally, we consider the effects of radia- 
tive losses on the spectrum in order to examine the spectral 
cutoffs produced, which are vital for comparison with obser- 
vations. It is clear that the complexity of such an algorithm 
would be significant and would require the use of current 
HFC methods. We now present an algorithm and associ- 
ated implementation designed to meet these requirements 
with the intention of proving the feasibility of this approach 
and providing a basis for further research. 



2 EQUATION, DISCRETISATION AND 
NUMERICAL SCHEME 

2.1 A relativistic transport equation 

Wc begin with the relativistic transport equation with units 
chosen such that c — 1. Initially we consider isotropic pitch 
angle diffusion in the co-moving (fluid) frame although this 
can be relaxed later with only minor adjustments to the 
model and the code. We will also assume that the hydro- 
dynamical structure of the shock is steady, though a time 
dependent model could be incorporated. The injection term 
Q is not explicitly defined here. In this analysis, any change 
to the injection term will be absorbed into Q. The stan- 
dard relativistic transport equation is expressed in a mixed 
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coordinate system as 



dx 



df 



+ ^(p)|((i-.^)|)-r(« + M)-,^ 



dp 



p'L{p)f) + Q. 



(1) 



Here the particle momentum p and direction angle /i are 
measured in the co-moving (fluid) frame, while time t and 

position X (and hence the background flow velocity u{x)) are 
measured in the shock rest frame. In general the pitch angle 
diffusion coeSicient D depends on the particle momentum. 
The energy loss rate L represents radiative losses and we 
will consider synchrotron radiation in this discussion. Hence 
g represents the angular behaviour of the synchrotron losses. 
Transforming the variable fi to the shock rest frame using 



u' — u(x) I a + u(x) 



(2) 



1 — n'u{x) " 1 + //m(x) ' 
we seek an equation for f{t,x,ij,',p). Firstly we can write 

, ^ /I — Ll'u fl'u — V? 

r 1 + MM = r r. — V + r- 

V 1 — M M 1 — A* w 

_ r _ 1 
~ r2(i -/i'ti) ~ r(i - /u'u) 

and it can be shown that 

(1-M)^ = (1-/. 

Prom these expressions one can derive the relativistic trans- 
port equation with //' measured in the shock rest frame, with 
the appropriate transformation applied to g to give g: 

d£ _^ ,df ^ ^2du ii'in' - u) d£ 
dt dx dx (1 — fj,'u) dp 



(3) 



(4) 



+ D(J>)T\l-u^l' 



dn' 



9/ 
dn' 



+ .9(M')r(l- V)^^ {pL{p)f)+Q (5) 

Let y = log(p/po)- For synchrotron losses we then have 
L{y) = Lople^^ (6) 



g{^J!) = 1 - /ttyi„id = 1 - ( r 

\ 1 Li 



= - . (7) 



Recovering factors of c where appropriate, and setting p = 
w/c, we can make the substitution / = /e"''^ to take advan- 
tage of the expected power law spectrum with index s ~ 4. 
Hence / 



dt dx (l-/i'/9) dx \dy \ 
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d_ 
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dl 
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r(l - fi'P) dy 



(8) 



We now consider the simplified case of a momentum in- 
dependent diffusion coefficient where D{y) = Dq ^ Q which 



is also constant across the shock. Do has units of inverse 
time, and we can divide equation 8 by Dq to find 

df ,d! _ -. ii'iiji' -p) dp (dl_.7 
dr di (1-At'^) di \dy ^ 

+ r3(i-,,'f|-((i-,'^)^ 

+ T^(^#-^f + ^ 
r(l - tx'l3) Do dy 



where 
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(10) 
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We are free to choose any momentum as the reference 
momentum po- We could choose the non-relativistic cutoff 
momentum associated with the synchrotron losses given by 



Po=P 



31/0 K 

Aw 



1 1 

— + — 

Wl W2 



(12) 



where Am = wi — M2 and the subscripts 1 and 2 repre- 
sent quantities measured upstream and downstream respec- 
tively. An expression describing the relationship between the 
pitch angle diffusion coefficient D and the spatial diffusion 
coefficient k, which is the basis of the diffusion approxi- 
mation, can be found in Kirk & Duffy (1999). Hence in 
the non-relativistic limit with D = Z)o(l — p^) we expect 
K = c^/{6Do) so that 



Po = 



Lo 



2DoA/3 V/3i ^2 



(13) 



where p = u/c. Therefore if we set the unitless parameter 



2/0 



■log 
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(14) 



(15) 



then wc have the non-dinicnsionalised equation 

9f I ,df ^ 2 f^'{l^' -P) dp /a/ _ - 
dr^^d^ (l-l^'p) d^ydy \ 

+ r3(i-,,r|((i-.'^)| 

, (1-m'^) Jy-yp) ^ , o 
+ r(l-M'/3)' dy^"^- 

Given an expression for the normalised flow speed ^(a;) and 
the injection term Q the solution of this equation can be 

scaled to find a specific solution for known values of Lo 
and Do. Furthermore, for P{x) <C 1, we expect the (non- 
dimensional) momentum cutoff to occur at j/ = 0. 

A very similar transformation can be applied for mo- 
mentum dependent diffusion. This will be the subject of fu- 
ture work. 



2.2 Notation 

For e£ise of use we will revert to the symbols x and t for the 
non-dimensional space and time variables and the transfor- 
mations outlined above are to be understood in their defi- 
nition. We will also dispense with the primes and bars on 
the quantities /i' and /. Therefore it should henceforth be 
understood that all quantities except for p are measured in 
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the shock rest frame, p being measured in the comoving fluid 
frame. We also note that / now describes deviations of the 
distribution function from a power law of index s. 

2.3 Discretisation 

The discretisation of the domain in jj, and y that we have 

chosen is very straightforward. In y a uniformly spaced grid 
of 0(1000) points is used. In /i we also choose a uniform 
grid with the spacing chosen to result in a few tens of grid 
points. The only subtlety involved is that we choose to place 
the first and last grid points a distance of A/x/2 inside the 
geometric end points fi — ±1. This choice greatly simplifies 
the boundary conditions for the chosen numerical scheme, 
as described in section 2.5. 

Regarding the x and t domains, two factors must now 
be accounted for: 

(i) The spatial domain must be greater than the difi^usion 
length scale at all (most notably high) momenta. 

(ii) The temporal domain must be greater than the dy- 
namical and acceleration time-scales. 

In addition, if we wish to compare our results to those ob- 
tained for discontinuous shocks, we should choose the shock 
width (and hence the spatial grid spacing) to be less than 
the diffusion length-scale at all (most notably low) momenta. 
Of course the method is not limited to such cases and in 
fact this constraint increases the computational complexity 
of the problem significantly. Let £{y) denote the diffusion 
length-scale at a given momentum and w the shock width. 
For an explicit finite difference scheme we then have 



cAt ~ Aa; < W < i{ymin) ^ £{ymax) <. Ctf. 



(16) 



These requirements can more easily be met through the use 
of a non-uniform grid spacing in x. 

One of the first tests of any code investigating DSA is 
the production of a power law with the appropriate spectral 
index. For a strong, non-relativistic shock (compression ratio 
of 4), with a momentum independent diffusion coefficient we 
expect a spectral index of —4. It has been shown (Kruells & 
Achterberg 1994) that a finite difference scheme modeling 
a continuous velocity profile will only produce the correct 
power law if 

(i) the advection length-scale per time step is less than 
the shock width and 

(ii) the diffusion length-scale per time step is greater than 
the shock width. 

The first of these conditions is trivially met for explicit finite 
difference schemes. The second condition requires considera- 
tion when choosing the ratio C of the shock width to the grid 
spacing at the shock. For a fixed shock width we can proceed 
to increase C by trial and error until the spectral index con- 
verges. A useful "rule of thumb" can be found by combining 
the above inequalities resulting in the requirement that 

2Atl3 
CAx 

where wq is the shock width in units of the diffusion length- 
scale. Hence as a starting point wo can choose the num- 
ber of points resolving the shock and then adjust the shock 
width to meet this condition. The situation is less clear for 



Wo < 



(17) 



relativistic flow where the diffusion length-scale is not well 
defined. 



2.4 Numerical scheme 

Finite difference methods are one of the most commonly 
used numerical approaches to solving partial differential 
equations. After years of development, a wide variety of 
schemes now exist to match the wide spectrum of prob- 
lems they are used to tackle. In general, the application of 
finite difference methods to a given equation requires care- 
ful selection of an appropriate approximation scheme from 
a relatively short list of commonly used ones. It is often 
instructive to start with the simplest, most crude schemes 
before deciding which of the more complex and accurate 
schemes could yield the most efficient results. We use the 
method of lines here, meaning that we will replace all deriva- 
tives w.r.t, {x,y,ii} with finite difference approximations 
and solve the resulting equation for ^ to advance the so- 
lution in time. This solution can be approximated by tak- 
ing discrete steps forward in time using a standard O.D.E. 
method. We have implemented the ubiquitous fourth order 
Runge-Kutta method, known for its reliability, and a third 
order Adams-Bashforth method. 

The simplest approach to the diffusion term is the stan- 
dard centred difference approximation to the second deriva- 
tive which is second order accurate in A/i. The source terms, 
containing no derivatives, do not require a difference approx- 
imation and are trivially implemented in the scheme. The 
only remaining terms are advective, representing the uni- 
form motion, shock acceleration and radiative losses (decel- 
eration) of particles. 

Because the MPI parallelisation has been implemented 
in the y dimension, we use the Lax-Wcndroff scheme in y 
to minimize communication costs while giving second order 
accuracy. A third order scheme is used in x to provide appro- 
priate accuracy for the steep gradients near the shock with- 
out incurring excessive computational cost. It behaves well 
on the non-uniform grid and while high order schemes can 
be unstable (especially where large gradients are present), 
this one seems well suited to our particular problem. The 
scheme for a non-uniform grid is given in the appendix. 

In summary, the numerical approximations used for 
each derivative in equation 15 are therefore the following 

^ : The third order Adams-Bashforth method (explicit, 
multi-step), or fourth order Runge-Kutta (multi-stage). 

1^ : The third order, upwind-biased scheme (see ap- 
pendix). 

■ The Lax-Wendroff scheme (second order). 



dy 



(jl — fj,^)^^ : Central differencing (second order). 



2.5 Boundary conditions 

The final requirement for the implementation of this scheme 
is a set of boundary conditions. Far downstream of the shock 
we expect the spatial gradient of the distribution function 
to be close to zero. Furthermore, any errors incurred there 
should have a minimal effect on the solution near the shock. 
The upstream boundary is more critical at the shock due to 
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the bulk flow direction. The distribution function decays ex- 
ponentially as we move upstream from the shock and, far up- 
stream, tends to zero. Thus the boundary condition should 
approximate either an exponential decay or could be set to 
zero if sufficiently distant from the shock. 

Suitable boundary conditions in pitch angle (cosine) 
are perhaps less obvious. In early testing a zero gradient 
condition was used, accurate to first order, at fj, = ±1. 
This did produce reasonable results, comparing well to 
the expected distributions except at the boundary points. 
While a second order accurate, zero gradient condition pro- 
duced slightly better results a more elegant solution was 
subsequently found. By choosing the grid points to lie 
at {-1 + A^i/2, -1 -I- 3Aij,/2, . . . , 1 - Aij,/2} in conjunc- 
tion with the standard central difference approximation to 
^^(1 — the values of the distribution function at 

jj, = ±1 conveniently vanish from the calculation. Thus no 
boundary condition is in fact required. Of course a boundary 
condition can be chosen, a posteriori, and fitted to the data 
if more detail near the end points is ever required. Alterna- 
tively we could extrapolate a suitable fitting function (sec 
section 4 for an example) to evaluate the solution very close 
to /u = ±1. The extrapolation distance can be reduced with 
higher resolution given sufficient computational resources. 

Finally, we must choose boundary conditions in momen- 
tum (magnitude). In the presence of losses it is clear that the 
solution will decay to zero in the limit of large momentum. 
Since the shape of the solution in the cutoff region is not 
generally known it is difficult to make any further assertions 
without reducing their applicability to special cases. Thus 
the upper boundary condition is chosen to lie well above 
the cutoff where a function value of zero is enforced. Ex- 
cluding points close to the shock, the lower boundary con- 
dition in momentum is somewhat less important because 
loss processes dominate any minimal acceleration and cause 
particles to flow outwards across the boundary. Indeed, test- 
ing has shown the lower boundary condition in momentum 
away from the shock to have negligible effect. The situation 
is somewhat complicated close to the shock where (for most 
pitch angles) acceleration can overcome losses leading to an 
inflow of particles across the boundary. This inflow of parti- 
cles is effectively an injection process, representing the point 
at which particles enter the domain. Since the detail of the 
injection process is not generally known, we have chosen to 
implement a constant source of particles at the shock with a 
sink (zero boundary condition) elsewhere on the lower mo- 
mentum boundary (which are of little significance to the 
solution) . 

3 DOMAIN OF INTEREST AND 

COMPUTATIONAL PARALLELISM 

The first task in the implementation of a finite difference 
scheme is to establish the domain and resolution of the nu- 
merical grid. In the x direction, the domain must enclose 
the region of plasma that can significantly affect the dis- 
tribution at the shock. Hence the downstream boundary is 
chosen to lie at Gl'iiy) where G ^ 1 and the subscript 2 de- 
notes a downstream quantity. The diffusion length-scale in 
the non-relativistic case is given by ^2(y) — K.{y) /u2 where 
K is the spatial diffusion coefficient and U2 the downstream 



fiow speed. Beyond the boundary position, particles are suf- 
ficiently unlikely to diffuse back to the shock that they can 
be neglected. In the upstream region the solution function 
is expected to decay exponentially with increasing distance 
from the shock. We therefore set our upstream boundary at 
a distance such that the particle phase space density has de- 
cayed to a negligible amount. A numerical value for the ap- 
proximate exponent associated with this decay can be found 
from the useful expressions given in Dempsey & Kirk (2008). 

The y domain can I)c chosen based on the estimated 
value of the momentum cutoff po, which implies a cutoff 
near y = 0. Based on previous work such as Dempsey & 
Duffy (2007), we can anticipate that the cutoff should take 
place over 2 or perhaps 3 decades in momentum. Hence, in 
order to examine the cutoff and power law region, we could 
choose the domain y £ [—10,3]. Of course this estimate can 
be revised if it becomes apparent that the necessary features 
are not captured. If we inject particles in a narrow energy 
range, it is also necessary to leave some room at the low 
energy end of the spectrum for the oscillations described in 
Kirk & Schneider (1989) to decay. 

The domain in /x must obviously span the range (—1, 1). 
The resolution must be sufficient to resolve any features in 
the distribution, but this scale is not immediately obvious 
a priori. We can, of course, run the code and decrease A/j 
as necessary for adequate resolution. Even in the case of 
high velocity flow, it appears that only a few tens of points 
are necessary for sufficient angular resolution. Of course this 
may not be the case if n is measured in a different frame. 

In order to attain more reasonable run times, paral- 
lelisation using MPI is necessary. The grid has largest di- 
mension in the y direction and hence we implement MPI 
parallclisation in that dimension. We use the Stokes cluster, 
run by the Irish Centre for High End Computing (ICHEC), 
which consists of 320 hex-core SMP nodes with ConnectX in- 
finiband interconnects. Non-blocking communications were 
employed, as well as parallel data output, using the MVA- 
PICH MPI implementation. 

Improved scalability can be attained on such SMP clus- 
ters through the use of a hybrid OpenMP-MPI implemen- 
tation. This is achieved by reducing the number of MPI 
processes in operation (which reduces the relative size of 
the communication halo and hence the comrrmnication over- 
head) and running numerous OpenMP threads per MPI pro- 
cess on each SMP node. Thus running N threads per pro- 
cess we reduce the relative size of the halo by a factor of A'". 
The total factor by which parallclisation overhead is reduced 
depends on the OpenMP overhead incurred, but for many 
problems it can be reduced by a factor of order N. 



4 RESULTS 

The primary test used in the development of this method 
was the production of a power law with appropriate spectral 
index. If we temporarily ignore the underlying physics that 
produces a particular compression ratio, we can investigate 
the spectral indices produced at a non-relativistic shock. We 
should find that the spectral index is given by the expression 
s = 3a/{<7 — 1) where a is the compression ratio. For test- 
ing purposes we chose an upstream fiow speed of ui = 0.1c 
with an isotropic, momentum independent diffusion coeffi- 
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Figure 1. Spectral Index vs Compression Ratio for a shock with 
upstream flow speed 0.1c. The line represents the non-relativistic 
analytical result in the diff'usion approximation. The crosses mark 
the numerical results tested. 
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Figure 2. The steady state distribution in momentum near a 
non-rclativistic shock of compression ratio 4 integrated over pitch 
angle (a) at the shock and (b) ~ 10 diffusion lengths downstream. 
Momentum independent, isotropic pitch angle diffusion and syn- 
chrotron losses used for comparison with Heavens & Meisen- 
heimer (1987). 

cient. No radiative losses were included for this test. The 
results arc shown in Figure 1 where we can see the method 
accurately reproduces the expected results. 

Figure 2 compares the results of our method to ana- 
lytical work by Heavens & Meisenheimer (1987) in the case 
of a non-relativistic shock with momentum independent dif- 
fusion. While the analytical results were obtained using the 
diffusion approximation, we can see that the results are quite 
accurate and it is difficult to differentiate between the curves 
in each pair above the injection energy. In this case the com- 
pression ratio of the shock was set at 4 resulting in a spec- 
tral index of 4.001 between the injection and cutoff energies. 
Hence the relative error in the spectral index produced by 
this method is less than 0.1% for the grid resolution used. 
The synchrotron cutoff region shows the strong agreement 
between our numerical results and the analytical ones of 



Heavens & Meisenheimer (1987). There is a slight discrep- 
ancy between the numerical and analytical results down- 
stream, the numerical cutoff being slightly sharper. This 
may be the result of the increasing spatial grid spacing. The 
results shown use a magnetic field proportional to the flow 
velocity and a diffusion coefficient that is constant in space. 

Near the injection energy the method also produces the 
oscillatory behaviour elegantly explained in Kirk & Schnei- 
der (1989), corresponding to the various generations of par- 
ticles rccrossing the shock after initial monocncrgctic injec- 
tion. Because the energy gain is pitch angle dependent for 
each interaction with the shock, the peaks in the spectrum 
spread out and settle down to a smooth power law. This in- 
jection region is not our primary interest in this study and 
hence we have not included the necessary mathematical de- 
tails to represent particles with speeds significantly less than 
c. Therefore we will not analyse this oscillatory behaviour 
any further, other than to note that some dampened oscil- 
lations are present at low energies as expected. 

For rclativistic shocks it is difficult to find purely an- 
alytical results for testing purposes. We turn to the eigen- 
function method for the most accurate results available that 
can provide useful tests for our code. Once again, our first 
test is a check on the spectral indices produced where radia- 
tive losses are unimportant. We have employed the Jiittner- 
Synge equation of state in order to calculate the compres- 
sion ratio. A similar scenario was considered by Kirk et al. 
(2000), where the effects of magnetisation were also taken 
into account in the compression ratio. Using Dexter, a 
data extraction applet provided by the Astrophysics Data 
System (ADS), reasonably accurate numerical values can 
be obtained for comparison with our results. Figure 3 com- 
pares the values extracted with those calculated using our 
method. It is clear that the values obtained for the spec- 
tral index are consistent, following this relatively complex 
curve quite effectively. Of course a very precise comparison 
would require more accurate knowledge of the results from 
the eigenfunction method and significant computational re- 
sources for higher resolution in our numerical results. 

Regarding strongly relativistic shocks, we note that eis 
the Lorentz factor of shock increases the rates of diffusion 
and acceleration increase rapidly. This introduces stiffness 
in the equation, the timcstcp being severely restricted by 
the factor of in the diffusion term. Ideally a new numer- 
ical scheme should be constructed to deal with this prob- 
lem. A split method using an implicit scheme or "super- 
timestepping" for diffusion could potentially be used. De- 
pending on the thickness of the shock profile, the factor 
in the acceleration term might also lead to stiffness, requir- 
ing additional splitting. The present implementation does 
not employ such methods, but works well for Lorentz fac- 
tors up to 10 or more. 

We would also like to verify that the method produces 
the expected behaviour at thick shocks. In order to do 
this we can compare the spectral index produced using our 
method with those produced by the eigenfunction method 
for thick shocks. Schneider & Kirk (1989) investigated this 
effect by assuming a power law solution without radiative 
losses in the steady state. Once again, we have used the 
Dexter data extraction applet to compare results from that 
paper with our numerical work in Figure 4. A shock of speed 
0.9c and compression ratio 2.43 is considered for a range of 
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Figure 3. Spectral index at a weaJcly magnetised thin shoclc com- 
pared witii results from Kirk et al. (2000) for a range of shock 
speeds. 
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Figure 4. Spectral index at a shock of speed 0.9c and compres- 
sion ratio 2.43 compared with results from Schneider & Kirk 
(1989) for a range of shock widths. Non-dimensional units are 
shown, as outlined in section 2.1. 



shock widths w. While the precision of the comparison is 
limited by the accuracy of the data extraction, we can see 
that the results agree quite well even for these expedient, low 
resolution tests. We should mention that our implementa- 
tion uses a error function expression to represent a sigmoid 
velocity transition at the shock instead of the hyperbolic 
tangent expression used by Schneider & Kirk (1989). Quali- 
tatively the two are very similar. Only a slight adjustment in 
the units was required to generate the plot shown in Figure 
4, and so we did not bother to change our shock profile for 
this comparison. 

Another test can be performed examining the angular 
distribution of particles at the shock. It has been shown that 
the first eigenfunction in the expansions used by Kirk et al. 
(2000) provides a reasonable approximation at relativistic 
shocks. Indeed, a good deal of further analysis can be per- 
formed using a single eigenfunction approximation as shown 
in Dempsey & Kirk (2008). With pitch angle measured in 
the shock frame and momentum (magnitude) measured in 
the fluid frame Kirk et al. (2000) give 



fs oc exp 



1 + H 



while Dempsey & Kirk (2008) give 



fs oc exp Our. ^ oc exp a„(l - pi] 



(18) 



(19) 



3 



0.6 



0.4 







New Method X 
Kirk et al. 
Dempsey & Kirk 


/Vx 
/V/ 
















% 

Xx 
Xx 

TNx^ 



Figure 5. Angular distribution at shock for ui = 0.9c 

where Uu is the solution of the transcendental equation 
exp(2au) = (a„(,5 - 1) — l)/(a„(/? + 1) — 1) (easily approx- 
imated numerically). Clearly the two eigenfunctions will be 
similar if Uu — (/3i — Using a numerical root find- 

ing algorithm it is straightforward to compare the two and 
show that the fractional difference falls below 10% near 
Pi = 0.5c, decreasing rapidly for faster shocks. Indeed, Kirk 
et al. (2000) state that their first eigenfunction (in isolation) 
only provides a good approximation above 0.5c whereas the 
eigenfunction used by Dempsey & Kirk (2008) yields accu- 
rate results at lower velocities. 

Figure 5 shows the angular dependence of the distribu- 
tion function at a shock of upstream flow speed 0.9c. It is 
clear that the two eigenfunction solutions are nearly iden- 
tical in this case. We can sec that our numerical method 
produces a very similar angular spectrum even from a low- 
resolution grid {Afj, ~ 0.143 in this case). Figure 6 shows the 
angular distribution at a shock of upstream flow speed 0.1c. 
Here we can clearly see that the eigenfunction of Dempsey & 
Kirk (2008) retains its accuracy at low velocities and agrees 
well with the numerical results. We expect the solution to 
become isotropic in the fluid rest frame at large distances 
downstream of the shock and this is also confirmed in the 
numerical results. Note that the end points of the numerical 
domain in lie within half of the grid spacing of the physical 
extrema n = ±1. Hence, for our numerical scheme, the van- 
ishing factor of 1 — /tt^ makes it unnecessary to evaluate any 
derivatives at = ±1, thus removing the need for boundary 
conditions. Note that we are not imposing the condition that 
there arc no particles at /i = ±1. If a boundary condition 
were necessary for a different scheme, a reflecting condition 
would be appropriate. 

Finally we show an example of a particle spectrum pro- 
duced by our method to demonstrate its use. Figure 7 shows 
the spectrum of high energy particles produced at a shock of 
speed 0.995c (F ~ 10) in the presence of synchrotron losses. 
The synchrotron loss rate used in the upstream region differs 
from that used downstream by the square of the shock com- 
pression factor, which is approximately 3. The diffusion co- 
efficient is constant in both space and momentum, rescaled 
to 1. The shock velocity transition is based on the error 
function on a length-scale of 10~^ in the units described in 
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Figure 6. Angular distribution at shock for ui = 0.1c 
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Figure 7. Particle spectrum for an ultrarelativistic shock of 
width 10~^ (see 2.1 for units). 



section 2.1, that being about 8 times longer than the de- 
cay length-scale of the solution upstream. Hence we are not 
modelling a "thin" shock here. The spectral index s in the 
power-law region (well below the cutoff) is approximately 
4.28. Also shown for comparison are two curves describing 
spectra of the form / oc p~"exp{—{p/pfit)'^) for ^ = 1,2 
where pfu was chosen to match the cutoff position. 



5 CONCLUSION 

We have presented a new approach for the solution of the 
particle transport equation about thick, hydrodynamically 
steady, plane relativistic shocks in one spatial dimension. 
This approach uses the fluid approximation to model evo- 
lution of the phase-space density function of the high en- 
ergy particles using finite difference approximations and 
the method of lines. Our objective is to combine and ex- 
tend the results obtained with the cigcnfunction method by 
various authors such as Kirk & Schneider (1987), Heavens 

6 Drury (1988), Kirk, Guthmann, Gallant & Achterberg 



(2000), Dempsey & Duffy (2007). Another significant ben- 
efit of the finite difference method of lines approach lies in 
the ease with which it could be coupled with a relativistic 
hydrodynamics code to model modified relativistic shocks 
(with time dependent structure). Our physical model ex- 
tends that of Kirk et al. (2000) with the addition of time 
dependence, a sigmoid shock velocity profile (for first order 
Fermi acceleration) and synchrotron losses. We have chosen 
a modified mixed coordinate system in which all quantities 
are measured in the shock rest frame except for p, the mag- 
nitude of the particle momentum, which is measured in the 
rest frame of the background (thermal) fluid flow. The angu- 
lar distributions produced using this coordinate system are 
quite smooth, which is rmmerically advantageous. 

The numerical scheme we have chosen is a fourth order 
Runge-Kutta integrator in conjunction with a second order 
finite diflFerence scheme in momentum and pitch angle and 
a third order scheme in space. A non-uniform grid is used 
in space to focus the computational effort about the shock. 
The scheme is implemented using a parallel, hybrid (MPI 
+ OpenMP) algorithm to take advantage of modern SMP 
cluster architectures, demonstrating good scalability up to 
many hundreds of CPUs. 

The method has been tested using analytical results in 
the non-relativistic limit. As expected, power law solutions 
are produced with spectral indices agreeing very well with 
standard theory for a range of compression ratios. At the 
synchrotron cutoff, the particle spectrum confirms the ana- 
lytical approximations of Heavens & Meisenheimer (1987), 
both at the shock and in the region around it. The angu- 
lar distribution is close to isotropy in non-relativistic cases, 
as we would expect, quantitatively reproducing analytical 
results. For relativistic shocks this method also produces 
power law solutions. The associated spectral indices can 
be compared with the numerical results obtained using the 
eigenfunction method, shown in Kirk et al. (2000). We have 
demonstrated that our method accurately reproduces these 
spectral indices from non-relativistic to ultra-rclativistic 
shocks. In all cases the angular distribution of particles was 
found to be in agreement with analytical approximations 
based on single eigenfunction solutions (Dempsey & Kirk 
2008). 

We have shown that standard CFD techniques can be 
used to find numerical solutions of the transport equation 
for high energy particles in the vicinity of relativistic shocks. 
A number of physical effects have been incorporated and it 
is possible, if not straight forward, to include additional ef- 
fects using these highly versatile finite difference methods. 
In addition to the results that can be obtained directly from 
this code, the methods employed would lend themselves par- 
ticularly well to combined hydrodynamical models in order 
to investigate acceleration at time dependent, modified rel- 
ativistic shocks. 
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APPENDIX A: FINITE DIFFERENCE 
APPROXIMATION SCHEMES 

The third order scheme wc have implemented can, for the 
constant velocity advection equation, be written 

U^T^ + (/l2 — hi)uT — hlh2 



fr=fr+uT 



{hohi + h^)h2 + hohl + + hi 



ft 



+ UT 



v^T^ + (/i2 — h\ — ho)uT + {—hi — ho)h2 
hohih2 + hohl 
t'^ + (h2 — 2hi — ho)uT 
{h\ + hohi)h2 



m 
Ji-l 



+ 



{{-2hi - ho)h2 + hi + hphi) 
{hi + hohi)h2 



fi 



M^r^ + {-2hi - ho)uT + {hi + hphi] 
hi + (2/ii + ho)hl + {hi + hohi)h2 



fr+i (Ai) 



where fto.i arc the nearest and next nearest grid spacings in 
the upwind direction and h2 is the nearest grid spacing in 
the other direction. 

Note that the Beam- Warming and third order schemes 
have an upwind bias so for m < apply the "reflection" 
transformation: hi —hi and fi±a fi^o,- 
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